# This file will is called by "MasterFile.R"
# remove(list=ls())
# setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# ------------------ Plot Specifications ---------------------
margins = c(3.5, 3.5, 1, 1)
plot_slippery_slope_two_economy = TRUE
plot_width = 4.2
plot_height = 3.8

# ----------------- Load the baseline calibrated parameters -------------
load(file = "baseline_calibration.Rdata")
# this loads the par vector for parameters

# calibration function 
source("Model scripts.R")  # This script solves the model and provide simulation functions 
calibration_results=calibration_moments(par)
equi_sol = calibration_results$equi_sol

# Then assign values for next processing
g_bar = par$g_bar; d_bar = par$d_bar; dt = equi_sol$dt;
AH = par$AH;  AL = par$AL
lambda=par$lambda; beta=par$beta; r = par$r; eta_H=par$eta_H; eta_L=par$eta_L
w_grid = equi_sol$w_grid
f_list = calibration_results$f_list
solve_individual_optimal_faster = f_list$solve_individual_optimal_faster; F_fun = f_list$F_fun; Delta_output_fun = f_list$Delta_output_fun; solve_equilibrium=f_list$solve_equilibrium; Simulate_model=f_list$Simulate_model; solve_other_moments=f_list$solve_other_moments; 
muK_fun=f_list$muK_fun; DeltaK_fun=f_list$DeltaK_fun; muw_fun=f_list$muw_fun; Deltaw_fun=f_list$Deltaw_fun; investment_cost_fun=f_list$investment_cost_fun;  l_dN_fun=f_list$l_dN_fun;  post_crisis_fraction_of_capital=f_list$post_crisis_fraction_of_capital; output_drop_and_gbar=f_list$output_drop_and_gbar;  
solve_welfare = f_list$solve_welfare
simulation_fun=f_list$simulation_fun
avg_zeta = par$avg_zeta; 
qL = equi_sol$qL; qH = equi_sol$qH;  

# Define capital drop and g_bar
Delta_capital_fun <- function(kappaH, kappaL, w){
  K_before = 1
  K_post = K_before * ( 1 + DeltaK_fun(kappaH, kappaL, w) )
  return( (K_post-K_before)/K_before  )
}
tol=10^(-6)
capital_drop_and_gbar <- function(g_bar, w_evaluation, qH = 2*beta+0.1, qL= beta+0.1, tol=10^(-6)){
  resultH = f_list$solve_individual_optimal_faster(qH, d_bar, g_bar)
  resultL = f_list$solve_individual_optimal_faster(qL, d_bar, g_bar)
  kappaH = integrate( function(zeta) f_list$F_fun(zeta+resultH$l_fun(zeta))*f_list$f_zeta_fun(zeta,"H"), lower=0, upper=zeta_upper, rel.tol=tol )$value 
  kappaL = integrate( function(zeta) f_list$F_fun(zeta+resultL$l_fun(zeta))*f_list$f_zeta_fun(zeta,"L"), lower=0, upper=zeta_upper, rel.tol=tol )$value 
  return( Delta_capital_fun(kappaH, kappaL, w_evaluation) )
}

# -------------------- Figure 8 with containing K obj: Slippery Slope in the Same Economy --------------------
# Government intervention and the avg w
g_bar_vec = c( seq(0, g_bar ,length.out=20), seq(g_bar+0.005, 0.25, length.out=20) )
index_baseline = which.min( abs(g_bar_vec-0.14)  )
# intervention and next intervention scale. 
govt_funding_next_crisis_vec = rep(0, length(g_bar_vec))
equi_sol_no_govt = solve_equilibrium(d_bar, g_bar=0)
equi_sol_no_govt = solve_other_moments(equi_sol_no_govt)
w0 = equi_sol_no_govt$w_avg
K0 = 1
# Choose a capital drop target that maps to no intervention at all. 
resultH = post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, 0, type="H")
resultL = post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, 0, type="L")
kappaH = resultH$kappa
kappaL = resultL$kappa
w = w0  + Deltaw_fun(kappaH, kappaL, w0) # immediately after the crisis
K = K0 * (1 + DeltaK_fun(kappaH, kappaL, w0) )
for(iter in 1:10/dt){
  # 10 years
  w = w + muw_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt
  K = K * (1 + muK_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt )
}  
capital_drop_target = capital_drop_and_gbar(0.096, w, equi_sol$qH, equi_sol$qL)  
# Think about the actual amount of intervention. 
govt_credit_provision_vec = g_bar_vec
govt_credit_provision_next_crisis = g_bar_vec
w_next_crisis_vec = g_bar_vec
K_next_crisis_vec = g_bar_vec
productivity_vec = rep(0,length(g_bar_vec))
productivity_next_crisis_vec = rep(0,length(g_bar_vec))
for(iter_g_bar in 1:length(g_bar_vec)){
  print( paste0("Progress:", iter_g_bar/length(g_bar_vec)*100, "%") )
  resultH = post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, g_bar_vec[iter_g_bar], type="H")
  resultL = post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, g_bar_vec[iter_g_bar], type="L")
  kappaH = resultH$kappa
  kappaL = resultL$kappa
  govt_credit_provision_vec[iter_g_bar] = w0*resultH$govt_sector_financing + (1-w0)*resultL$govt_sector_financing
  productivity_vec[iter_g_bar] = w0*AH + (1-w0)*AL
  w = w0  + Deltaw_fun(kappaH, kappaL, w0) # immediately after the crisis
  K = K0 * (1 + DeltaK_fun(kappaH, kappaL, w0) )
  for(iter in 1:10/dt){
    # 10 years
    w = w + muw_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt
    K = K * (1 + muK_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt )
  }  
  w_next_crisis_vec[iter_g_bar] = w
  K_next_crisis_vec[iter_g_bar] = K
  result = uniroot( function(g_bar) capital_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - capital_drop_target, interval = c(0, 2)  )
  govt_funding_next_crisis_vec[iter_g_bar] = result$root
  
  resultH= post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="H")
  resultL= post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="L")
  # Total credit provision
  govt_credit_provision_next_crisis[iter_g_bar] = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
  productivity_next_crisis_vec[iter_g_bar] = w*AH + (1-w)*AL
}
sensitivity = ( govt_funding_next_crisis_vec[index_baseline] - govt_funding_next_crisis_vec[1] ) / ( g_bar_vec[index_baseline] - g_bar_vec[1] )
pdf(paste0("../Figures/FigureA13a_same_economy_slippery_slope_w.pdf"), width=plot_width,height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( g_bar_vec, w_next_crisis_vec, xlab=TeX("current govt program size $\\bar{g}$"), ylab=TeX("next crisis (in 10 years) $\\omega$'") )
dev.off()
pdf(paste0("../Figures/FigureA13b_same_economy_slippery_slope_credit.pdf"), width=plot_width,height=plot_height, family = "serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( g_bar_vec, govt_funding_next_crisis_vec, xlab=TeX("current govt program size $\\bar{g}$"), ylab=TeX("next crisis (in 10 years) govt intervention $\\bar{g}'$") )
dev.off()



